# Purpose: Produce the county-level maps in the paper (e.g., Figure 2)

#####################################
## 1. Read in and prepare the data ##
#####################################

# See README for instructions on how to access this data
dat <- read_csv("../replication_regional_data/raw/map_measures/regional_measures.csv")

# Read in the shape files
shapes_in <- st_read("../replication_regional_data/raw/map_shapefiles/NUTS_RG_20M_2016_4326.shp",
                     layer="NUTS_RG_20M_2016_4326")

# Make nuts3 border data for Germany
nuts3_shapes <- shapes_in %>%
  filter(LEVL_CODE == 3) %>%
  filter(CNTR_CODE == "DE")

# Make nuts1 border data for Germany
nuts1_shapes <- shapes_in %>%
  filter(LEVL_CODE == 1) %>%
  filter(CNTR_CODE == "DE")

# Combine shapes and data
dat_map <- dat %>%
  filter(sample == "all") %>%
  select(nuts3, measure, sy_avg_resid, native_avg_resid) %>%
  left_join(nuts3_shapes, by=c("nuts3"="NUTS_ID")) %>%
  st_as_sf()

# Make the palette
pal5 <- (brewer.pal(n = 5, name = "YlGnBu"))

# Select the outcomes
outcomes <- c("n_frnd_nat_lcl", "n_lcl_50n_grps", "produ_any_de", "n_frnd_sy_lcl", "n_frnd_eu_lcl", "n_frnd_t5_lcl")


###########################
## 2. Make baseline maps ##
###########################

for(i in 1:length(outcomes)){
  
  curr_map <- dat_map %>%
    filter(measure == outcomes[i]) %>%
    mutate(sy_avg_resid_ntile = paste0(as.character((ntile(sy_avg_resid, 5) - 1) * 20), "-",
                                       as.character(ntile(sy_avg_resid, 5) * 20)))
  
  ggplot(curr_map) +
    geom_sf(data=nuts3_shapes, fill="gray", linewidth=0.02) +
    geom_sf(aes(fill = as.factor(sy_avg_resid_ntile)),  linewidth=0.02) +
    geom_sf(data=nuts1_shapes, fill=NA, linewidth=1) +
    theme_void() +
    labs(fill = "%-tile") +
    scale_fill_manual(values = pal5)
  
  ggsave(str_interp("../results_regional_analysis/nuts3_map__unshrunk__${outcomes[i]}.png"), last_plot(), width=5, height=5)
  
}


########################################
## 3. Make decomposition measure maps ##
########################################

dat_map <- dat_map %>%
  mutate(relative_friending = sy_avg_resid/native_avg_resid)

outcomes <- c("n_frnd_nat_lcl", "n_lcl_50n_grps")

for(i in 1:length(outcomes)){
  
  curr_map <- dat_map %>%
    filter(measure == outcomes[i]) %>%
    mutate(relative_avg_ntile = as.factor(paste0(as.character((ntile(relative_friending, 5) - 1) * 20), "-",
                                                 as.character(ntile(relative_friending, 5) * 20)))) %>%
    mutate(native_avg_ntile = as.factor(paste0(as.character((ntile(native_avg_resid, 5) - 1) * 20), "-",
                                               as.character(ntile(native_avg_resid, 5) * 20))))
  
  plot_measures <- c("relative_avg_ntile", "native_avg_ntile")
  
  for(j in 1:length(plot_measures)){
    
    ggplot(curr_map) +
      geom_sf(data=nuts3_shapes, fill="gray", linewidth=0.02) +
      geom_sf(aes_string(fill = plot_measures[j]),  linewidth=0.02) +
      geom_sf(data=nuts1_shapes, fill=NA, linewidth=1) +
      theme_void() +
      labs(fill = "%-tile") +
      scale_fill_manual(values = pal5)
    
    ggsave(str_interp("../results_regional_analysis/nuts3_map__${outcomes[i]}__${plot_measures[j]}.png"), last_plot(), width=5, height=5)
    
  }
}